Effects of imperfections for Shor's factorization algorithm 
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We study effects of imperfections induced by residual couplings between qubits on the accuracy of 
Shor's algorithm using numerical simulations of realistic quantum computations with up to 30 qubits. 
The factoring of numbers up to N = 943 show that the width of peaks, which frequencies allow 
to determine the factors, grow exponentially with the number of qubits. However, the algorithm 
remains operational up to a critical coupling strength e c which drops only polynomially with log 2 N. 
The numerical dependence of e c on log 2 N is explained by analytical estimates that allows to obtain 
the scaling for functionality of Shor's algorithm on realistic quantum computers with a large number 
of qubits. 
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I. INTRODUCTION 



It is quite clear that the main interest to the quantum 
computing, has been generated by Shor's factorization al- 
gorithm [l| which has exponential efficiency gain com- 
pared to any known classical algorithm. Indeed, Shor's 
algorithm allows to find the factors of a large number N 
with 0(ln N) quantum gates while all known classical 
algorithms require a number of operations that grows ex- 
ponentially with In TV (see e.g. review and relevant Refs. 
in 0). Without any doubt this result has a fundamental 
importance from a mathematical view point. However, 
its implementation in real life requires an understand- 
ing of effects of imperfections and errors unavoidably 
present in any real physical realization of the algorithm 
on a realistic quantum computer. Again, here a math- 
ematician can be satisfied by a mathematical statement 
that in quantum computations the errors grow not faster 
than quadratically with the number of performed quan- 
tum gates (see 0]) and thus the global accuracy of the 
algorithm is good enough if the norm of errors in each 
quantum gate is sufficiently small. However, a physicist 
generally would like to see more concrete and realistic es- 
timates of the algorithm accuracy. Unfortunately, direct 
experimental verification of the accuracy for a large num- 
ber of gates and qubits is not possible at present. Indeed, 
the most advanced quantum computation of Shor's algo- 
rithm has been done on a 7-qubit NMR-based quantum 
computer that allowed to factorize only a rather small 
number N = 15 (even if certain simplifications of the 
original algorithm have been used) [|[. 

Therefore, the only possibility remaining is the method 
of numerical simulations testing various types of realistic 
errors and imperfections. The first steps in this direction 
have been done in 0, H, @] • A number of interesting ef- 
fects of errors on the accuracy of Shor's algorithm has 
been found in these pioneering works but the factorized 
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number was still N = 15 and therefore it was not pos- 
sible to determine the accuracy scaling at large values 
of N. More recently, additional numerical studies have 
been performed to investigate the effects of finite accu- 
racy in quantum phase rotations of the quantum Fourier 
transform (QFT) algorithm used in Shor's factorization 
@, Hii dynamical phase errors in Shor's algorithm with 
N up to 33 @ and discrete qubit flip errors with N 
up to 247. In the latter case the QFT part of Shor's al- 
gorithm has been performed in a semiclassical way using 
the one qubit control trick (see e.g. [Til. Ill fill. fli| ) while- 
the modular multiplication has been been performed with 
up to 20 qubits including the workspace using the circuit 
described in [l5| . 

In this work we perform extensive numerical simula- 
tions investigating effects of imperfections on the accu- 
racy of Shor's algorithm factorizing numbers up to a max- 
imal value N = 943 using up to L = 30 qubits. We con- 
centrate our studies on the case of static imperfections 
which induce static one-qubit energy shifts and residual 
static couplings between qubits following the lines started 
in [l6l |. This type of imperfections is especially impor- 
tant since generally the errors produced in this case are 
accumulated coherently and lead to a more rapid drop 
of fidelity and accuracy of quantum computations com- 
pared to the cases of noisy unitary errors in qua ntum 
gates [13, [HI and dissipative decoherence [H, l20L HH . 
It is known that for the quantum algorithms simulating 
problems of quantum chaos periodic in time the Floquet 
eigenstates are exponentially sensitive to static imperfec- 
tions [22j . Due to that the study of their effects on the 
accuracy of Shor's algorithm becomes especially relevant 
since recently it has been shown that certain blocks of 
the Shor algorithm are characterized by the properties 
of quantum chaos [23| . Also it is important to note that 
Shor's algorithm is essentially based on a determination 
of a certain frequency of return. In some cases, like in the 
Grover algorithm, such a frequency can be exponentially 
sensitive to static couplings [24[ that makes the investi- 
gation of static imperfections effects in Shor's algorithm 
even more important. 

Thus, in the present work we present the first numcri- 
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cal studies of how the static imperfections affect the ac- 
curacy of Shor's algorithm. Our aim is to determine the 
parametric dependence of the accuracy on the imperfec- 
tion strength, number of qubits and number of gates. 
For this we use a simplified but generic model of im- 
perfections which can be applied to various implemen- 
tations of Shor's alg orithm discussed in the literature 
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The paper has the following structure: Section II gives 
a brief description of ideal Shor's algorithm, Section III 
describes the model of errors introduced by static imper- 
fections, the results of numerical studies are presented in 
Section IV and the discussion of the results is given in 
Section V. 



II. IDEAL REALIZATION OF SHOR'S 
ALGORITHM 

First we briefly describe the main structure of Shor's 
algorithm [l[ factorizing a large integer number N using 
parallelism of many-body quantum evolution. Following 
Shor we choose a random number x relatively prime to 
N and calculate its order r (also called period) defined as 
the minimal positive integer value such that 



1 mod N. 



(1) 



Once r is known there is a high probability to obtain 
two non-trivial factors of N by a classical computation 
in polynomial time (in the number of binary digits of TV). 
This procedure fails in rare cases [l[ and in such a case 
one has simply to chose a different value of x and restart 
again. 

The difficult task is to compute the order r and this 
task can be efficiently achieved by Shor's algorithm pro- 
vided we have a reliable quantum computer with a suf- 
ficient number of qubits at our disposal. This algorithm 
requires an i-qubit state composed of two quantum regis- 
ters which we will call the control register (with ni qubits) 
and the computational register (with n q = L — ni qubits). 
We associate to the basis states of both registers integer 
numbers by: 



10 = \a n -i)n-i ■ ■ ■ ■ ■ \ao)o 
where in binary representation 

I = a a + 2ai + . . . + 2™- 1 a„_i 
and n = ni for the control register or n = 



(2) 



(3) 



for the com- 



putational register. Here 



represents the jth qubit of 



the register and <Xj € {0, 1}. In order to factorize a num- 
ber N one needs to choose ni and n q such that 2™" > N 
and Q = 2 n ' > N 2 , therefore typically n ; sa 2n q . 
We first prepare the initial state 

l^) = |0)«,|l)n, (4) 



and then apply single qubit Hadamard gates to every 
qubit in the control register and get (dropping subscripts) 



(5) 



The principal idea of Shor's algorithm is the observa- 
tion that one can construct a combination of quantum 
gates, acting on both registers, that performs for all 
a = 0, . . . , Q — 1 simultaneously the operation : 



|a)|l) -► |a)|a; a mod N) 
which gives the state 



(0) 



(I) 



Then, after obtaining the state \1jj2), we apply the QFT 
[3] to the control register 

Q-1Q-1 

^ = nEE e l2 ™ c/Q \c)\x a mod N) (8) 

^ c=0 a=0 

and measure both arguments to get 



P(c, x ) = | (ip3\c}\x mod AO | 

2 



1 E 



ca/Q 



Q 



a:x" —x ' 



(9) 



where k = 0, . . . , r — 1 is arbitrary and the sum over a 
runs over all values such that x a = x mod N . There- 
fore a = rv + k where v = 0, . . . , — 1 and Mfc = 
\{Q — k — l)/r] + l and the evaluation of the sum yields: 



P(c,x k ) 



1 s\TL 2 {M k ircr/Q) 
Q 2 sin 2 (ncr/Q) 



(10) 



This function only depends weakly on the choice of k 
(since Q > N 2 and N > r > k such that Q 3> k and 
Mfc ~ Q/r > 1 is nearly constant in k) and as a function 
of c it has r equidistant strongly localized peaks of width 
unity, of height M 2 /Q 2 w 1/r 2 and located at mQ/r 
with m = 0, 1, . . . , r — 1. 

If the algorithm is run by an ideal quantum computer, 
then with a very high probability the outcome of a mea- 
surement will be given by an integer value of c which is 
very close to one of the peaks mQ/r. Thus, using a con- 
tinuous fraction expansion we can determine the rational 
number p/q closest to c/Q with a denominator smaller 
than N . Here the choice Q > N 2 ensures that there is 
at most one such number inside the peak and therefore 
p/q coincides with m/r. Furthermore the position num- 
ber m of the peak is quite random and if by chance m is 
relatively prime to r one obtains directly r — q and the 
algorithm succeeds. However, if m and r have a common 
divisor larger than unity we have q — r/ gcd(m, r) < r 



and the algorithm did not succeed. Therefore one has to 
check classically if the candidate u q" for r is indeed a so- 
lution of x q = 1( mod N). This fortunately can be done 
in a polynomial time. In case of failure the algorithm has 
to be repeated and even though the probability of suc- 
cess is not very high one obtains after a few (C(log log r)) 
measurements [l[ the correct value q = r. 

Practically it is more convenient to measure only the 
control register which provides c with the total probabil- 
ity: 

P{c) = Y J P{c^ k )~rP{c,x k ). (11) 

k 

We note that the dependence of M\. on k in Eq. (fTU|) is 
rather weak and therefore the above procedure to deter- 
mine the minimal period r remains the same. 

However, this description of the algorithm still lacks 
some precision how to implement the operation described 
in equation j6|). Suppose we are able to perform on the 
computational register the multiplication by x mod N: 

\y) -> U nlult (x)\y) = \(yx) mod TV) (12) 

by some unitary operator. Of course this operator can- 
not be unitary if we require this for all values y — 
0, . . . 2 nq — 1 simply because the classical application 
y — ► (xy) mod TV is not unique on this set (unless 
TV = 2" q which is of no interest). If we require that x 
and TV are relatively prime then this application is unique 
at least for y = 0, . . . N - 1 and for y = TV, . . . , 2 n - - 1 
we have to complete it in some unique way, for example 
by: y — > y if y > TV. Therefore wc define the quantum 
multiplication operator by x mod TV by: 

U^WW} - \ \y) f y = N,...,2 n * -1 

(13) 

The states \y) with y > TV are in principle not relevant 
for the ideal Shor algorithm because they are never pop- 
ulated in the perfect computation and the effect of the 
quantum gates on these states is rarely discussed in the 
literature [jj, la, l2a] • However, they are important to en- 
sure overall unitarity and they may be very well popu- 
lated if the quantum computation is subjected to errors 
or imperfections. Furthermore, we note that in the defini- 
tion (|13[) we could in principle replace the unit-operator 
acting on the non-relevant states by an arbitrary uni- 
tary operator (acting on a space of dimension 2 Tlq — TV) 
provided that we do not mix relevant (y < TV) and non- 
relevant states (y > TV). 

We now introduce the controlled multiplication op- 
erator Uq^ uU (x) acting on both registers (control and 
computational register) applying the simple multiplica- 
tion (fT3"|) on the computational register if and only if 
the j—th qubit of the control register is |1). Develop- 
ing a = Y^jLo a o 2-' with aj E {0, 1} we see that the 
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operation ([6]) can be performed by the unitary operator: 

ni—l 

t/ F mult(.T) = J] ^Cmult Hiod N) (14) 

3=0 

since 

ni—X ni — l 

*°=nKr= n * 2J ^ 

j=0 j=0,a j= l 

and where in the last equation every multiplication is 
taken modulo N. 



x a mod N 
I 1 




FIG. 1: Quantum circuit of Shor's algorithm on an ideal 
quantum computer with the quantum multiplication opera- 
tor Uj = f/ mu it {x 23 mod JV) as denned in (|13[) 

Fig. Q] shows the schematic quantum circuit of Shor's 
algorithm on an ideal quantum computer in terms of the 
quantum multiplication operator (fT3")) . To complete an 
explicit implementation one has to show that this op- 
erator can be realized in terms of elementary one or 
two qubit quantum-gates. We do not enter into de- 
tails here and mention as examples the important works 

0, [25[ that provided explicit implementations of the 
quantum multiplication by x mod N using 0{n^) ele- 
mentary gates. These implementations require also ad- 
ditional work space qubits which are initially |0...0) 
and must remain so after completion of this operator, 

1. e. the implementations must eventually provide code 
to reversibly "erase" the additional work space qubits. 
We assume that there are no errors inside this additional 
work space and no errors coupling it to the control and 
computational registers. In this way we may restrict our 
consideration only to L = n q + n\ qubits. 

III. SHOR'S ALGORITHM WITH STATIC 
IMPERFECTIONS 

We now turn to Shor's algorithm in the case of static 
imperfections generated by residual couplings be- 
tween qubits and energy level shifts. The effects of these 
imperfections and their numerical modeling have been 
considered in detail in [l8| on examples of quantum chaos 
algorithms (see also Refs. in [l~8| on other works). There 
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it has been shown that effects of static residual couplings 
can be modeled by an additional unitary rotation act- 



s'' mod N with static errors 



ing between two arbitrary gates : U s 



ASH 



Here SH 



represents the Hamiltonian due to the residual static cou- 
plings between qubits which provides a non-trivial evo- 
lution of the state stored in the quantum register even in 
absence of any quantum gate. In this approach the quan- 
tum gates are considered to be exactly ideal. In principle 
SH may couple all qubits in the control register, in the 
computational register and in the additional work space 
necessary for the concrete implementations of the quan- 
tum multiplication operator (|13p . However, in this work 
we use a simplified error model in which SH couples only 
the qubits in the computational register and therefore 
in Shor's algorithm the initial Hadamard gates or the fi- 
nal quantum Fourier transform are not affected by these 
errors. In principle the quantum Fourier transform is 
considered as relatively stable with respect to errors [H 
and the number of Hadamard gates ni is relatively small. 

Furthermore we do not consider a specific implemen- 
tation of the quantum multiplication operator (|13p . we 
only assume that it can be written as a product 



U 2 -U 1 



(16) 



where Uj , j = 1 , . . . the elementary quantum 

gates which constitute this operator and n m = 0(n^) is 
the number of these elementary gates. A specific choice 
of Uj depends on the classical variable x, and also on N, 
and since the value of x significantly affects the algorithm 
implementation we have a different set of gates Uj for 
each x (and N). 

Thus, in presence of static imperfections the quantum 
multiplication operator U mu \t(x) has the form: 



Umnit{x) = U„ 



ASH 



j j iSH j j iSH 

U2 ■ e ■ Ui ■ e 



(17) 



We now introduce an effective perturbation operator for 
the full multiplication operator by: 

U mult (x) = U mult (x)e* 6H ° !! ^ . (18) 

From Eq. (| 1 T[) we may determine SH c fi(x) as: 



iSH ea (x) _ A6H(n m -l) . 



JSH(l) , ASH 



with 



8H(j) = U 



j'-i 



C/f 1 SH Ux 



(19) 



(20) 



We mention that the precise relation between SH c g(x) 
and SH is not really important in our approach since we 
directly model SH c f[(x) in our numerical simulations and 
use the expression (fT8|) without entering into details of a 
particular implementation of U mu \t(x). We remind that 
in (|19[) . (|20|) the dependence of 5H c s(x) on x is given by 
the choice of elementary gates Uj which are changed with 
a change of x. A schematic quantum circuit of Shor's 
algorithm on a quantum computer with static imper- 
fections in the quantum multiplication operator (|13p is 
shown in Fig. [21 



ii>, 



u a - 



QFT 



-U x - 



FIG. 2: Quantum circuit of Shor's algorithm on a quantum 
computer with static imperfections in the quantum multipli- 
cation operator Uj = (7 mult [x 21 mod N) lfT3|) and the ef- 
fective static perturbation STLj = SH c g (x 21 mod N) [see 
Eqs. (|13|l . (HI}], where in numerical simulations STLj is given 
by Eq. (|21[) with random realizations of Si, Ji for practically 
each j (see text). 



As we already mentioned, the implementations of 
the quantum multiplication operator (|13p described in 
Refs. [5], [25[ require additional work space qubits which 
arc initially placed in the state |0...0) and are erased 
after the computation. The implementation of Ref. || 
contains a quantum code that erases the work space 
qubits correctly but only for the relevant states \y) with 
< y < N as initial states and not for the non-relevant 
states with y > N. For a perfect quantum computer this 
is of course not a problem, but when taking into account 
errors the non-relevant states may be populated and dif- 
ferent implementations, which are absolutely equivalent 
for the relevant states, may potentially behave quite dif- 
ferently with errors. Even if the particular implemen- 
tation ensures that a non-relevant state as initial state 
produces a properly erased work space register, the er- 
rors may still produce non-erased contributions. 

Actually the use of work space qubits implies that the 
notion of non-relevant states has to be enlarged, i.e. a 
combined state \y) IV'work) in the computational and work 
space register has to be considered as non-relevant if ei- 
ther y > N for iV'work) = 1 ... 0) or y arbitrary for 
l^work) -L |0...0). If Shor's algorithm is implemented 
on a perfect quantum computer without any imperfec- 
tions these non-relevant states are never populated. How- 
ever, errors and imperfections will populate these states 
and their role is potentially quite important in this con- 
text. In this work we do not want to enter into the details 
of the effects due to the work space qubits. So, we simply 
assume that our model of imperfection effects (fTT)) acts 
only in the computational register, or in other words the 
static imperfections do not couple computational qubits 
with work space qubits. However, even in this approx- 
imation we still keep track of the non-relevant states in 
the computational register (the states \y) with y > N). 

For numerical simulations of Shor's algorithm in pres- 
ence of imperfections we use a classical computer taking 
into account the control register (with up to 20 qubits) 
and the computational register (with up to 10 qubits) 
and up to 30 qubits in total. We do not implement the 
quantum multiplication operator in terms of elementary 
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gates but we directly implement the unitary operator as 
given in Eq. (|13jl . To model the static imperfections we 
used the multiplication operator with errors given by (fTF 
and with the effective perturbation operator given by: 



n„ — 1 



SH cS (x) = > 5,a. 



where are the Pauli operators acting on the ith qubit 
(of the computational register) and Sj, Jj are random 
coefficients, chosen differently for each value of x and 
distributed according to: 



H 

i=0 



0) 



2 V Ji^M+i ( 21 ) 



i=0 



5 h J, e [VSe, V3e] 



(22) 



We remind that even for static imperfections 8H c s(x) 
given by Eqs. (fT9| . ([20)) strongly depends on the actual 
value of x because this factor is hardcoded in realistic 
implementations by the choice of elementary gates Uj. 
According to Eq. ([H]) . we have to apply the (controlled 
version) of the multiplication operator for all values 



X G {x v mod N | j = 0, . . . , m — 1} 



(23) 



In our numerical simulations we have ensured by the 
proper choice of Si, Ji that SHc^x 2,3 mod N) is identi- 
cal to 5H e f[(x 2 mod N) if for j ^ I we have x 2,3 = x 2 
mod N, Otherwise, we have chosen different realizations 
of Si , Ji for each value of x 2J mod N assuming that the 
x-dependence of the hard coded implementation is suffi- 
ciently complex to render SH c g(x) uncorrelated for dif- 
ferent values of x. This introduces some kind of slight 
correlation that takes into account the static property of 
the imperfections. However, we have also checked that 
neglecting these correlations (choosing each time a dif- 
ferent realization of SH c ff(x) even if the same x- value ap- 
pears again) does not affect significantly our numerical 
results discussed below. We also note that in potential 
applications (for "real" quantum computers) with larger 
values of TV and n; , n q these kind of correlations will be- 
come less important. So, in majority of cases for each j 
we have <5-ff e ff(j) with independent random realizations 
of 5i, Ji in Eq. (|2"Tj) distributed as in (|A"| . 

In principle the unitary operator e lSHee{ ( x > is quite ran- 
dom due to Eqs. (fl9|) . ([20]) and should directly couple 
many qubits in the computational register. Our model 
(|2ip for the imperfections is quite convenient for numer- 
ical computations and is similar to the model used in 
[l7L [IH but with a difference that in those works it is 
the elementary residual Hamiltonian SH which is chosen 
in this way [see Eq. I|17[)]. Therefore we should expect 
that SH e s(x) has a more complicated structure than ([21) . 
However, choosing Sj and Jj of comparable size we are 
well in the quantum chaos regime [lg, EH, HH and there- 
fore the model ([2~Tj) describes well the effects of static 
imperfections. It should be also noted that the quantum 
gates of the algorithm introduce additional strong mixing 
between all qubits even if they are not directly coupled 
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FIG. 3: Probability P(c) (Eqs. (dTTJ) at the final stage of 
Shor's algorithm for e = (top) and e = 0.1 (bottom) for 
values N = 323, n q = 9, L = 27, x = 2, r = 72. 



by SH e s(x) (see [22( for details). We also checked that 
a change of SH e g(x) from the form of Eq. ([21]) to the 
case when all qubits are coupled by residual interactions 
does not affect significantly the results of numerical sim- 
ulations. This is in agreement with the results obtained 
in [TH, [HI. Another advantage of a choice of 8H e s{x) 
in the form ([2"Tjl is the local structure of couplings be- 
tween qubits that corresponds to a physical reality. It is 
also important to note that when we have n g = rii gates 
as in Fig. [2l then the effective strength of e is effectively 
renormalized as e — * e^/TT^ since static errors in each real- 
ization of 8TLj (see Fig. [2|) are independent and random. 
We leave the question about possible strong correlations 
between STlj due to a specific implementation of the al- 
gorithm for future studies. 

The above consideration assumes that sufficient ran- 
domization of static imperfections takes place along the 
path of a specific quantum circuit for the modular multi- 
plication. In this case we may assume that the effective 
Hamiltonian SH c g(x) in the propagator contains differ- 
ent random couplings between qubits for each value of 
x (see Eqs. (|19HAI0 . However, it is possible that the 
errors remain well correlated along the path of this cir- 
cuit and in this case it is more appropriate to consider 
that SH c g(x) does not depend on x and remain the same 
along the whole Shor's algorithm. In our numeral studies 
we mainly concentrate on the first possibility ("generic 
imperfection model") but in order to have the complete 
picture of the effects of static imperfections we also con- 
sidered the second case with 5H c r(x) remaining constant 
along the full circuit ("correlated imperfection model"). 
According to our previous discussion the important prop- 
erty of both models is that the errors appear only via the 
positions of the propagator e l8H " a ^ between the mod- 
ular multiplications in the full circuit of the algorithm. 
Hence the specific implementation of the modular mul- 
tiplication circuit does not affect the random properties 
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of inter-qubit couplings in SH c g(x). Therefore, once the 
parametric dependence on the imperfection strength e, 
number of qubits n q and number of gates n g is estab- 
lished through numerical simulations, we can apply these 
results to arbitrary implementations currently discussed 
in the literature. 

In some sense, the model of static errors considered 
here can be viewed as a kind of generic static error model. 
It shows sufficiently rich and generic effects of errors and 
due to its certain simplicity allows to make numerical 
simulations with factorization of larger N values com- 
pared to previous numerical studies [f| 0, E] ■ This 
allowed us to determine the accuracy dependence on the 
parameters and to obtain the scaling law for a large num- 
ber of qubits. This required to perform extensive numer- 
ical simulations with up to 30 qubits which became pos- 
sible because we neglected the errors in the work space 
qubits. However, as soon as we obtain the parametric de- 
pendence of the algorithm accuracy we may reincorporate 
the effect of imperfections in the work space by modifying 
the effective qubit number in the computational register. 
We also neglected the static imperfections in the control 
register since the number of gates in the QFT (operating 
in the control register) is much smaller than the number 
of gates in the main part of Shor's algorithm. However, in 
the case of the correlated imperfection model, we verified 
that the introduction of couplings in the control register 
does not modify the established parametric dependence 
on the number qubits. We emphasize that our numer- 
ical calculations keep the exact quantum entanglement 
for the whole quantum evolution with up to 30 qubits. 
We note that a further increase of the factorized num- 
ber N can be achieved by replacing the control register 
by one qubit combined with appropriate measurements 
of this qubit and a semiclassical implementation of the 
QFT [lfl, EL EE El- However, this approach simu- 
lates the quantum measurement process in the algorithm 
and does not give a direct access to the full probability 
distribution in the quantum register which is substan- 
tially used in our studies. 

We present obtained numerical results in the next Sec- 
tion. 



IV. NUMERICAL RESULTS 

The effects of static imperfections in Shor's algorithm 
are studied numerically following the approach described 
in the previous Section: a wave vector of size 2 L is 
propagated numerically according to the quantum cir- 
cuits shown in Figs. 1 1121 all quantum gates are assumed 
to be exact, the imperfections, induced by residual cou- 
plings between qubits in the computational register, are 
encountered by the propagators exp(iSHj) appearing n; 
times in the circuit as it is described in Fig [2] We fac- 
torize numbers TV up to N f=a 1000. This means that 
we simulate numerically a quantum computer with up to 
30 qubits, 10 computational qubits and 20 control qubits 



(we assume ideal evolution in the workspace) . The list of 
factorized numbers N used for numerical simulations is 
given in Table I. We try to consider mainly most difficult 
cases when N has only two factors and their values are 
more or less comparable. 



TABLE I: Values for the data presented in Figs. HE] and [TO] for 
the generic imperfection model. Only the values with symbols 
are plotted in Figs. 18191 
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In Fig. [3] we show a typical example of the probability 
distribution P(c) of Eq. (jTTJ) for the ideal case e = 
(top) and for e = 0.1 (bottom). It can be seen that the 
imperfections significantly reduce the amplitudes of the 
main r peaks and lead to appearance of new small peaks 
in new positions. 

Since the success of the algorithm depends essentially 
on a probability of hitting r— peaks in the process of mea- 
surement then the most direct way to study this prob- 
ability is by clashing all the peaks into one, or in other 
words, adding them all together by taking c modulus 
s where s is the nearest integer value of the ratio Q/r 
and thus reducing all probabilities inside one cell with s 
states. In this way we obtain a new distribution of global 
search probability W(c): 

r-l 

W(c) = p (i c + s + JQM mod s) (24) 

3=0 

where now c = —s/2,...,s/2 — 1 (the difference of c 
for P and W is clear from the context) and s « Q/r 
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FIG. 4: (color online) The global probability distribution 
W(c), as denned in Eq. I)24J) . averaged over 10 realizations 
of random static imperfections, for different values of cou- 
pling strength e : (a) e = 0.025, (b) e = 0.05, (c) e = 0.1, (d) 
e = 0.2. The fast oscillating green (gray) lower curve shows 
the theoretical probability G(c) at e = (Eq. (f25)) . The solid 
(black) curve in (a) is the actual probability at e = ob- 
tained numerically. The red (dark gray) curves show W{c) 
obtained numerically at given values of e > 0. Here, as in 
Fig. [3] N = 323, n q = 9, L = 27, x = 2 and r = 72. 



is the distance between peaks. For the ideal algorithm 
this global probability W(c) has one peak at c = that 
stresses the important property of Shor's algorithm: it 
is not important what peak from the main chain of r 
peaks is selected by measurement, it is important to know 
its exact position modulus s that allows to determine 
r value and then to find the factors of N by classical 
computations. The global probability W(c) is distributed 
over states with c=— s/2, ... ^ s/2 — 1 and is normalized 
to unity in this interval. 

In Fig. S] we show a typical example of the global 
probability W(c) variation with the increase of coupling 
strength e. The distribution W(c) for the ideal algo- 
rithm at e = is well described by the envelope func- 
tion Wq(c) = (sin(7rc)/(7rc)) 2 of the distribution G{c) 



discussed in [31| (see also (fTO 



G(c) 



sin(7rc) 



sin(7rcr/Q) 



(25) 



Shor's algorithm is successful if the probability at c = 
is significant (comparable to 1). This is indeed the case 
for small values of e (Fig. Hh,b) . In these cases the main 
probability is concentrated near c = 0. There are new 
peaks appearing at very large values of c but they have 
rather small total probability. With a further growth of 
e the number of such peaks and their probability grow 
(Fig. 3J:), the amplitude of the central peak at c = 
drops and above certain e the distribution W(c) becomes 
practically flat (Fig. [4ji) that signifies the complete de- 
struction of the algorithm. A pictorial view of variation 
of W(c) with e is shown in Fig. O 
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FIG. 5: (color online) Quantum melting of Shor's algorithm 
induced by imperfections: color density plot of the global 
search probability W(c) as a function of coupling strength 
between qubits e for N = 323, n q = 9, L — 27, x = 2, r — 72 
(W(c) is averaged over 20 realizations). 

In order to study the effects of static imperfections on 
the algorithm accuracy in a more quantitative way it is 
convenient to use the inverse participation ratio (IPR) 



(26) 



which gives a number of effectively populated states in 
the distribution W(c). This quantity is extensively used 
to characterize the prop erties of many-body quantum 
states (see e.g. [HI, l32|). Another convenient charac- 
teristics is the width of the distribution defined as 



An= Ef(c)(c-(c)) ! 



(27) 



The dependence of these quantities on the perturbation 
strength e is shown in Figs. [S] and [7] for the typical case 
TV = 323. The value of £ is practically constant up to a 
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value e c « 0.1 after which it starts to grow abruptly. On a 
contrary, the width An grows starting from small values 
of e. At large e the saturation of growth takes place due 
to a finite number of states inside the distribution W(c). 




c 



FIG. 6: Averaged IPR £ given by Eq. (l26jl as a function of 
e for N = 323, n q — 9, L — 27, x = 2, the inset shows the 
dependence on small e in log-scale, average is done over the 
number of realizations given in Table I. 




e 



FIG. 7: Same as in Fig. El but for An given by Eq. p7)l. 

The dependence of the IPR £ on e for different N is 
shown in Fig. [8l The data clearly show that the depen- 
dence becomes more and more sharp with the increase of 
N. For An we see a strong increase with N but there is 
no such sharp behavior (see Fig. [9]) . We attribute such 
a difference to the fact that even small e gives far tran- 
sitions with exponentially large c ~ s sa Q/r oc 2 nq . Due 
to that the second moment of the probability distribution 
grows exponentially with the number of qubits. The nu- 
merical data on dependence of An on N at a small fixed 
e indeed show the exponential growth with An r* AeN 
with a numerical constant Aw 14 (Fig. [9] bottom panel). 
A similar behavior has been seen for quantum chaos al- 
gorithms [33l . |34| . The mechanism of this exponential 
growth is the following [H, [34| : the gates with imperfec- 
tions transfer a probability W e ~ e 2 n q from the search 
state at c w to about n q peaks (see Fig. [5]) distributed 
in the interval of size s ~ N . Here, n q comes from the 
norm of the Hamiltonian (|21[) with n q qubits with local 



couplings. There are n g — ni m 2n q such transitions 
W e during the whole algorithm computation. Thus, we 
obtain the second moment of the distribution W(c): 

(An) 2 w a 2 e 2 n q N 2 , (28) 

where according to numerical data of Fig. [5] (bottom) 
the numerical coefficient a f=a Aj \Jn q ' sa 4.5 is close to 
the one obtained in [H, [HI- Of course, the fluctuations 
in Fig. [9] (bottom) are rather large. We think that the 
main origin of these fluctuations is related to the arith- 
metic properties of x, r and N. Indeed, r varies signifi- 
cantly with x and N (see Table I) that clearly affects the 
transition probability induced by imperfections (3o| . In 
spite of these fluctuations the global exponential growth 
of An with n q is seen rather clearly. Such an exponen- 
tial sensitivity of (An) 2 on N is not very pleasant for 
the algorithm accuracy, but in principle this behavior is 
not so dangerous. Indeed, the total probability to have 
exponentially large values of c is very small and doing 
a few measurements and making a majority "vote" will 
eliminate such extreme values of An. 




0.1 0.2 0.3 
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FIG. 8: (color online) Dependence of averaged IPR £ on 
strength of imperfections e for different values of N, curves 
with symbols from the list of Table I (from top curves with the 
largest iV = 943 to bottom curves with the smallest N = 14); 
number of disorder realizations used for averaging is given in 
the Table. For the large values of N we show typical statisti- 
cal error bars, for small N the error bars are comparable with 
the symbol size and we do not show them. Top and bottom 
panels show £ in normal and logarithmic scale respectively. 

Therefore, more crucial is the behavior of £ since above 
a certain critical value e c the probability W{c) spreads 
over very many levels and the algorithm stops to work. 
Indeed, it is known that static imperfections can lead to 
a complete derealization, for example, in the case of a 
quantum algorithm simulating the Anderson localization 
in three dimensions 13611. 
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FIG. 9: (color online) Top panel: dependence of averaged 
An on e for different TV with the same symbols from Table I 
as in Fig. [8] Bottom panel shows the dependence of An on 
TV in log-log scale for e = 0.05, the straight line shows the 
dependence An = AeN with A w 14. 
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FIG. 10: Dependence of e c , obtained from the criterion 
£(e c ) = 10£(e = 0), on log 2 TV in log-log scale. The numerical 
data are shown for the generic imperfection model (squares, 
data from Table I), the correlated imperfection model with 
qubit couplings in the computational register (full circles) and 
the correlated imperfection model with all qubits in the con- 
trol and computational registers coupled by interactions (open 
circles), see text for model description. The straight lines 
show the fit e c = S/(log 2 TV)' 5 in the interval 4 < log 2 TV < 10 
with B = 0.98 ± 0.16, /3 = 1.04 ± 0.094 for squares (top line), 
B = 2.06 ± 0.42, /3 = 1.6 ± 0.11 for full circles (middle line) 
and B = 0.33 ±0.05, /3 = 1.57 ±0.09 for open circles (bottom 
line) . 



To determine the derealization border for Shor's al- 
gorithm and the dependence of e c on TV in the generic 
imperfection model we use a numerical criterion £(e c ) = 



10£(e = 0). Indeed, an increase by a factor 10 is suffi- 
ciently large to obtain the transition border in e. The 
dependence of e c on TV is shown in Fig. [TO] (squares and 
top line). From the theoretical view point the errors 
are accumulated randomly so that the probability Wt 
transferred from c = to all other states grows propor- 
tionally to the number of gates n g with errors and thus 
W t ~ W t n g ~ e 2 n q n g ~ e 2 n g . We expect that above 
the border Wt ~ 1 the probability becomes delocalizcd 
over exponentially many states and the algorithm is de- 
stroyed. This gives the quantum chaos border 



; (TV) = 73/log 2 (TV)« V2B/^T- 



(29) 



where B is a numerical constant. For our generic im- 
perfection model we have in the second equality n g w 
2n q but in the case when the errors related with the 
workspace qubits are taken into account we have n g ~ 
n q ~ (log 2 TV) 3 . The numerical data for e c are pre- 
sented in Fig. [TO] The fit of the dependence in the form 
e c = B/(\og 2 N) !3 in the interval 4 < log 2 TV < 10 gives 
B = 0.98 ± 0.16, (3 = 1.04 ± 0.094. Thus, the numerical 
data confirm the theoretical estimate (|29|) with B w 1. 
The border e c drops polynomially with log 2 TV since the 
whole Shor algorithm is performed in a polynomial num- 
ber of gates n g ~ (log 2 TV) 3 . In this respect the situation 
is different from the case of the Grover algorithm with im- 
perfections considered in [2l| where the number of gates 
grows exponentially with n q . The exponential sensitivity 
of Floquet eigenstates to static imperfections in quantum 
chaos algorithms [22j also corresponds to a different sit- 
uation since in a sense eigenstate corresponds to a very 
long time scale where the number of gates becomes ex- 
ponentially large. 

In the above consideration for the generic imperfection 
model, we assumed that the quantum circuit effectively 
modifies the couplings between qubits in the propagator 
from one gate to another. Another limiting case corre- 
sponds to the correlated imperfection model, where these 
couplings remain unchanged from gate to gate (see Sec- 
tion Hn]). For this particular model we also performed 
extensive numerical simulations considering two cases a) 
the interactions exist only between qubits in the compu- 
tational register (see Fig. [TQl full circles, middle line) and 
b) the interactions exist between all qubits in the control 
and computation registers (see Fig. [TQl open circles, bot- 
tom line). For the numerical study of these two cases we 
used the same quantities as those described above for 
the generic imperfection model. We do not reproduce all 
data here but only show the cumulative final dependence 
for the quantum chaos border e c defined by the same re- 
lation £(e c ) = 10£(e = 0). The fit of the numerical data 
in the form e c = 73/(log 2 TV)^ gives the same exponent 
^ w 1.6 for both cases of the correlated imperfection 
model with the numerical factors as in Fig. [TO] Nat- 
urally B becomes smaller when all qubits are coupled. 
The value of (5 is definitely larger as compared to the 
generic imperfection model (where f3 w 1). This can be 
understood on the following physical grounds: the errors 
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accumulate coherently along n g gates so that the transi- 
tion probability from the target state to all other states 
is W t ~ W e n^ ~ e 2 n q n 2 ~ e 2 n q . The quantum chaos 
border is given by the condition Wt ~ 1 that gives: 

e c (A0 = 5/ log 2 (A0 3 / 2 « ; (30) 

since we always chose n g ~ 2n 9 . The theoretical expo- 
nent /? = 1.5 is in good agreement with the numerical fit 
[3 = 1.6±0.1. We also clearly see that the fact of coupling 
all qubits does not affect the parametric dependence of 
the chaos border on log 2 (iV) and gives only a change of 
the numerical prefactor B. It is important to note that 
the quantum chaos border is lower for the correlated im- 
perfection model. 



between qubits and a coupling strength implementing a 
two-qubit gate. According to the above estimate in a 
quantum computer this ratio should be kept as small as 
Jrcs/Jg < £c ~ 2 x 10~ 6 to have a possibility to beat 
a modern classical computer in the RSA-factorization. 
Such a restriction rises serious requirements to experi- 
mental implementations of quantum computers, but it's 
possible to hope that future technological progress will 
make this possible. Finally we note that we do not con- 
sider here quantum error corrections (see [38[ and Refs. 
therein) which may improve the situation but on a price 
of significant increase of the total number of qubits re- 
quired for computations. 

This work was supported in part by the EC IST-FET 
project EuroSQIP. For numerical simulations we used the 
codes of Quantware Library [39j . 



V. CONCLUSION 

We performed extensive numerical simulations of 
Shor's algorithm factorizing numbers up to N = 943 
on a quantum computer with up to 30 qubits in pres- 
ence of residual static couplings between qubits. Our 
studies show that the width An of r-peaks, which po- 
sitions are essential for determination of factors of N, 
grow exponentially with N (see Eq. (|28|) ). However, the 
use of majority vote with few measurements allows to 
eliminate the rare events which contribute to this ex- 
ponential growth. In fact the algorithm remains op- 
erational up to the critical coupling strength e c which 
drops polynomially with log 2 N (see Eq. (|29|) ). Since 
with the work space qubits the total number of gates in 
Shor's algorithm is n g ~ (log 2 Nf the relation (J29J) gives 
e c ~ l/(log 2 A) 2 . In this estimate, based on Eq. (|2"9"|) 
with n q ~ log 2 N and n g ~ (log 2 -/V) 3 , we assume the va- 
lidity of the generic imperfection model where couplings 
fluctuate from gate to gate. Another limit corresponds 
to the case of correlated imperfection model where cou- 
plings remain fixed for all gates. In this case the re- 
lation (j30|) gives e c ~ l/(log 2 N) 7 / 2 . A presence of fi- 
nite correlation length 1 < 

figcor — in the number of 
gates n g will give interpolation between these two limit- 
ing cases with e c ~ l/[(log 2 N) 2 \jn gcor \. At present, the 
latest RSA challenge number factored is RSA-640 with 
log 2 A r = 640 [13 ■ Thus, assuming a more optimistic 
case of the generic imperfection model, a quantum com- 
puter which factors this number should have a dimen- 
sionless coupling strength e < e c ~ 2 x 10~ 6 . The value 
of e can be interpreted as e w J reB St, where J ros is a 
strength of residual couplings and St « 1/ J g is a time 
duration of two-qubit gate which is related to a typical 
value of coupling J g between two qubits which imple- 
ments this gate. As a result, we obtain that e ~ J Tes / J g 
has the meaning of the ratio between a residual coupling 



APPENDIX A: ERRATUM 

In the published version [Phys. Rev. A 75, 052311 
(2007)] as well as in the original arXive version there is 
an error in the numerical implementation of the static 
imperfection model (Eq. (21) in the paper) which re- 
sulted in an effective reduction of the disorder strength 
e by approximately a factor of 2 for the cases of generic 
and correlated imperfection models. However, there was 
no error for the correlated imperfection model with all 
qubits in the control and computational register coupled 
by interactions. After correction of the error the depen- 
dence of the inverse participation ratio (IPR) £ on e (sec 
Fig. [TT] (a) ) remains qualitatively the same as in Fig. 6 
of the paper. Nevertheless, the values of e c fluctuate 
strongly depending on the arithmetic properties of x, r 
and N. These fluctuations remain quite strong even af- 
ter increasing the number of data points (see Table II and 
Fig. [TT1 (b) as compared to Table I and Fig. 10). The al- 
gebraic fit e c = £/(log 2 Nf gives In 5 = 0.068 ± 0.105, 
j3 = 1.420 ± 0.054 for the generic imperfection model, 
InB = 0.70 ± 0.19, (3 = 1.897 ± 0.097 for the cor- 
related imperfection model, and ln£> = —1.22 ± 0.13, 
/3 = 1.523 ± 0.068 for the correlated imperfection model 
with all qubits coupled. We attribute the deviations of 
the exponents (3 from the theoretical values to finite size 
effects in number of qubits and the significant arithmetic 
fluctuations of e c as a function of r, x and N . 

In addition after publication we realized there were two 
typos in Eqs. (22) and (26) which should read 



8 U Ji€[-V3e,V3e] ; £= £ \W{c)\ 2 




This did not affect any of the calculations. 
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TABLE II: Values for the data presented in Fig. [TT] O) with 
tc^ for the generic imperfection model (□), e^ 2 ' for the cor- 
related imperfection model(») and ei^ for the correlated im- 
perfection model with all qubits coupled (o). Nr is the ap- 
proximate number of random realizations for each case. 
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